#!/bin/bash
# multi gvcf to vcf

#GATK dir
gatk=/picb/lilab/tools/gatk-4.2.6.1/gatk

#Reference dir
reference=/picb/lilab5/share_data/ref/hg38/
bundle=/home/lilabguest7/lilab5_cxf/ref/gatk_bundle

#shell options
working_dir=$1
sample_list=$2
merge_name=$3

samples=$(cat ${sample_list})

if [ ! -d ${working_dir}/vcf ]
then mkdir -p ${working_dir}/vcf
fi

if [ ! -d ${working_dir}/info ]
then mkdir -p ${working_dir}/info
fi 

if [ ! -d ${working_dir}/vcf/merge_vcf ]
then mkdir -p ${working_dir}/vcf/merge_vcf
fi 

for i in {1..22} X Y ; do
    sample_gvcfs=""
    for sample in $samples ; do
        sample_gvcfs=${sample_gvcfs}"-V ${working_dir}/gvcf/${sample}.HC.g.vcf.gz "
    done
    time $gatk CombineGVCFs \
        -L $i \
        -R $reference/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
        ${sample_gvcfs} \
        -O ${working_dir}/vcf/${merge_name}.HC.${i}.g.vcf.gz \
        > ${working_dir}/info/4.gvcf_to_vcf.info && \
    time $gatk GenotypeGVCFs \
        -L $i \
        -R ${reference}/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
        -V ${working_dir}/vcf/${merge_name}.HC.${i}.g.vcf.gz \
        -O ${working_dir}/vcf/${merge_name}.HC.${i}.vcf.gz \
        >> ${working_dir}/info/4.gvcf_to_vcf.info
done && wait
merge_vcfs=""
for i in {1..22} X Y; do
    merge_vcfs=${merge_vcfs}"-I ${working_dir}/vcf/${merge_name}.HC.${i}.vcf.gz "
done && time $gatk MergeVcfs ${merge_vcfs} -O ${working_dir}/vcf/merge_vcf/${merge_name}.HC.vcf.gz

if [ ! -d ${working_dir}/vcf/vqsr ]
then mkdir -p ${working_dir}/vcf/vqsr
fi 

#VQSR
time $gatk VariantRecalibrator \
    -R ${reference}/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -V ${working_dir}/vcf/merge_vcf/${merge_name}.HC.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 ${bundle}/hapmap_3.3.hg38.ensembl.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 ${bundle}/1000G_omni2.5.hg38.ensembl.vcf.gz \
    --resource:1000G,known=false,training=true,truth=false,prior=10.0 ${bundle}/1000G_phase1.snps.high_confidence.hg38.ensembl.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 ${bundle}/Homo_sapiens_assembly38.dbsnp138.ensembl.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
    -tranche 100.0 -tranche 99.8 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
    -mode SNP \
    --tranches-file ${working_dir}/vcf/vqsr/${merge_name}.HC.snps.tranches \
    --rscript-file ${working_dir}/vcf/vqsr/${merge_name}.HC.snps.plots.R \
    -O ${working_dir}/vcf/vqsr/${merge_name}.HC.snps.recal && \

time $gatk ApplyVQSR \
    -R ${reference}/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -V ${working_dir}/vcf/merge_vcf/${merge_name}.HC.vcf.gz \
    --truth-sensitivity-filter-level 99.8 \
    --tranches-file ${working_dir}/vcf/vqsr/${merge_name}.HC.snps.tranches \
    --recal-file ${working_dir}/vcf/vqsr/${merge_name}.HC.snps.recal \
    -mode SNP \
    -O ${working_dir}/vqsr/${merge_name}.HC.snps.vcf.gz \
    >> ${working_dir}/info/4.gvcf_to_vcf.info \

time $gatk VariantRecalibrator \
    -R ${reference}/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -V ${working_dir}/vqsr/${merge_name}.HC.snps.vcf.gz \
    --resource:mills,known=false,training=true,truth=true,prior=12.0 ${bundle}/Mills_and_1000G_gold_standard.indels.hg38.ensembl.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 ${bundle}/Homo_sapiens_assembly38.known_indels.ensembl.vcf.gz \
    -an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
    -mode INDEL \
    -tranche 100.0 -tranche 99.95 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
    --tranches-file ${working_dir}/vqsr/${merge_name}.HC.indels.tranches \
    --rscript-file ${working_dir}/vqsr/${merge_name}.HC.indels.plots.R \
    -O ${working_dir}/vqsr/${merge_name}.HC.indels.recal && \

time $gatk ApplyVQSR \
    -R ${reference}/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -V ${working_dir}/vqsr/${merge_name}.HC.snps.vcf.gz \
    --truth-sensitivity-filter-level 99.95 \
    --tranches-file ${working_dir}/vqsr/${merge_name}.HC.indels.tranches \
    --recal-file ${working_dir}/vqsr/${merge_name}.HC.indels.recal \
    -mode INDEL \
    -O ${working_dir}/vqsr/${merge_name}.HC.VQSR.vcf \
    >> ${working_dir}/info/4.gvcf_to_vcf.info

bgzip ${working_dir}/vqsr/${merge_name}.HC.VQSR.vcf && \
mv ${working_dir}/vqsr/${merge_name}.HC.VQSR.vcf.gz ${working_dir}/vqsr/${merge_name}.HC.VQSR.vcf.bgz && \
tabix ${working_dir}/vqsr/${merge_name}.HC.VQSR.vcf.bgz
